#! usr/bin/perl -w
use strict;

# print results on two output files
sub output
{
	my($ref_list_in,$ref_list_out)=@_;
	my @list_in=@$ref_list_in;
	my @list_out=@$ref_list_out;
		
	open(OUTPUT1,">SNP_in.txt") or die "cannot create outputfile";
	open(OUTPUT2,">SNP_out.txt") or die "cannot create outputfile";

	foreach my $item1 (@list_in)
	{
		print OUTPUT1 "$item1\n";
	}		
	
	foreach my $item2 (@list_out)
        {
                print OUTPUT2 "$item2\n";
        }

	close OUTPUT1;
	close OUTPUT2;

	print "\nsee independent SNPs in the output file \"SNP_in.txt\" and others in \"SNP_out.txt\" \n\n";

}

sub inLD
{	
	my $snpTest=shift;
	my($ref_data,$ref_list_in)=@_;
	my %data=%$ref_data;
	my @list_in=@$ref_list_in;	
	my $item="";
		
	#print "le contenu de list_in : "."@list_in";
	
	print "\n--------------- sub inLD ----------------\n";
	print "\nle snp a tester : $snpTest \n";	
	if((scalar(@list_in)) ne 0)
	{
		print "le contenu de list_in : "."@list_in";
		print "\n les snps dans list_in : ";
		
		foreach $item (@list_in)
		{
			print $item;
			while ( my($k,$v) = each(%data))
                        {
				print "\n dans le while item : $item, k : $k, v : $v \n";
                               	if( (($k eq $item)and($v eq $snpTest)) or (($k eq $snpTest)and($v eq $item)) )
                               	{
				       return "true";
                               	}
                               	else
                               	{
                                       return "false";
                              	}
			}
		}
		#for(my $i=0; $i<=(scalar(@list_in)); ++$i)
        	#{
		#	print $list_in[$i]." ".$i;
        	#	while ( my($k,$v) = each(%data))
                #	{
                		#if( (($k eq $list_in[$i])and($v eq $snpTest)) or (($k eq $snpTest)and($v eq $list_in[$i])) )
                 #       	{
                  #      		return "true";
                   #     	}
                        #	else
                        #	{
                         #       	return "false";
                        #	}
       			# }
#		}	
	print "----------- fin sub inLD ------------";
	}
	else 
	{
		return "false";
		print "----------- fin sub inLD ------------";
	}
}

# return two lists : list_in and list_out(SNPs'list to delete)
sub selection
{
	my($ref_list_in,$ref_data)=@_;
	my @list_in = @$ref_list_in;
	my %data = %$ref_data;
	my @list_LD=();	

	foreach my $key (sort keys(%data))
	{
		my $val = $data{$key};
               	my @key_in=grep(/^$key$/ig,@list_in);#search key in list_in (list of independent SNP)
		my @key_LD=grep(/^$key$/ig,@list_LD);#search key in list_out		

		print "\n\n ####### Key $key has value $val #########\n";
		
		if((@key_in eq 0) and (@key_LD eq 0))#key doesn't exists in list_in and list_out
		{
			my $IsinLD = &inLD($key,\%data,\@list_in); #SNP in LD with another in list_in ?
			print "\n$key en LD avec un snp dans list_in ? $IsinLD \n";
			push(@list_in,$key) if($IsinLD eq "false");
		} 
	}
	
	my @all_snp=(keys%data,values%data);
	my @uniq_all_snp = keys %{{ map { $_ => 1 } @all_snp }};
	my @list_out=();
	
	# create list_out
	foreach my $item (@uniq_all_snp)
        {
                push(@list_out,$item) if((grep(/^$item$/ig,@list_in)) eq 0 );
        }
		
	return(\@list_in,\@list_out);
}

sub load_data
{
	# parameters
	my ($filename,$t)=@_;

	my $firstLine=1;
	my @list_in=();
	my %data=();

	open(LDFILE, $filename) or die "cannot open $filename";
	while(my $line=<LDFILE>)
	{
		chomp($line);
		$line=~ s/ +/ /gi;#separate rows with a whitespace only
		my @fields=split(/ /,$line);

		if($firstLine!=1)#not header
		{		
			# collect fields of interest
	                my $R2=$fields[(scalar@fields)-1];# field R2
        	        my $SNPid1=$fields[3];# SNP1
                	my $SNPid2=$fields[6];# SNP2
			
			# r2 < threshold : SNP independent so put in list_in
			if ($R2<$t)
			{
				my @match1=grep(/^$SNPid1$/ig,@list_in);#search SNPid1 in list_in 
				my @match2=grep(/^$SNPid2$/ig,@list_in);#search SNPid2 in list_in 
				
				push(@list_in,$SNPid1) if(@match1 eq 0);
				push(@list_in,$SNPid2) if(@match2 eq 0);			
			}
			# else r2 > threshold : put in table to analyse
			else
			{
				$data{$SNPid1}=$SNPid2;
			}
		}

		$firstLine=0;

	}
	close LDFILE;

	# SNP selection
        my($ref_new_list_in,$ref_list_out)=&selection(\@list_in,\%data);	
	my @new_list_in=@$ref_new_list_in;
	my @list_out=@$ref_list_out;
	
	return(\@new_list_in,\@list_out)
}


sub main{

if (scalar@ARGV != 2)
{
	print "\nExpected 2 parameters. \$ perl snpIndep.pl <LDfilename> <threshold>\n";
}
else
{	
	# parameters
	my $LDfilename=$ARGV[0];
	my $threshold=$ARGV[1];
	
	# load data and SNP selection
	my($ref_list_in,$ref_list_out)=&load_data($LDfilename,$threshold);
	my @list_in=@$ref_list_in;
	my @list_out=@$ref_list_out;

	# output
	&output(\@list_in,\@list_out);
}


exit;
}

main();
